Self-organized pattern formation and noise-induced control from particle computation 
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We propose a new non-equilibrium model for spatial pattern formation on the basis of local 
information transfer. Unlike standard models of pattern formation it is not based on the Turing in- 
stability. Information is transmitted through the system via particle-like excitations whose collective 
dynamics result in pattern formation and control. Here, a simple problem of domain formation is 
addressed by this model in an implementation as stochastic cellular automata. One observes stable 
pattern formation, even in the presence of noise and cell flow. Noise stabilizes the system through 
the production of quasi-particles that control the position of the domain boundary. Self-organized 
boundaries become sharp for large system with fluctuations vanishing with a power of the system 
size. Pattern proportions are scale-independent with system size. Pattern formation is stable over 
large parameter ranges with two phase transitions occurring at vanishing noise and increased cell 
flow. 



An astonishing property of the development of mul- 
ticellular organisms is its extreme error tolerance and 
robustness against perturbations A key mechanism 
that coordinates structure formation during development 
is the self-organization of spatiotemporal patterns of 
gene activity [J- While detailed dynamical models of 
the involved gene regulation processes are not within 
close reach (mostly due to lack of genomic details) , phe- 
nomenological models of developmental processes have 
been studied for quite a while. One standard model for 
pattern formation in development is diffusion-driven pat- 
tern formation exploiting the Turing instability 3] . This 
principle has been applied to modeling biological organ- 
isms 0,0 and is able to account for a number of observed 
features of developmental processes, e.g., in the fresh wa- 
ter polyp Hydra |(|. 

However, diffusion based models, dating from the pre- 
genomic era, are by definition limited. When applying 
this type of model to an organism as well known as Hy- 
dra, for example, several unresolved problems persist. 
E.g., morphogen molecules postulated by the model still 
have not been identified. Further, parameter fine-tuning 
is needed, including a non-trivial hierarchy (separation 
by orders of magnitude) of diffusion constants |jj • Most 
importantly, several experiments point at specific devel- 
opmental features that are not easily captured by these 
models as, e.g., the reorganisation of the body pattern 
from a fully random cell assembly Q and the extreme 
sharpness and vertical regulation of expression bound- 
aries 0, including scale- invariant proportions. While a 
chemical gradient along the body axis could in principle 
provide position information, the precision of "gradient 
readout" is low at typical (i.e. low) morphogen concentra- 
tions This indicates that the information processing 
in the gene regulation machinery during development is 
not accurately captured in the diffusion-based picture. 
Indeed, recent experiments show that a large number of 
regulatory genes are involved in cell differentiation and 
pattern formation [Til ] . It may therefore be time to think 
about alternatives to diffusion-driven pattern formation, 



and in particular explore the possibilities of information- 
flow-driven pattern formation. The aim of this paper is 
to contribute to this discussion. 

Let us first recapitulate a simple developmental prob- 
lem, then define a simple stochastic cellular automata 
model that solves this pattern formation task. It is first 
demonstrated that this system performs de novo pattern 
formation, independent of initial conditions. Then sta- 
bility of pattern formation is studied in the presence of 
noise and cell flow. 

Position-dependent gene activation is a frequently ob- 
served mechanism in animal development. One example 
is the fresh water polyp Hydra, which has three distinct 
body regions - a head with mouth and tentacles, a body 
column and a foot region. The positions of these regions 
are accurately regulated along the body axis. In addi- 
tion, new cells continuously move from the central body 
region along the body axis towards the top and bottom, 
and differentiate into the respective cell types according 
to their position on the head- foot axis. This cell flow 
requires considerable robustness of the regulatory pro- 
cesses. The two remarkable features of this regulation 
that we focus on are the scale-independent position reg- 
ulation and the "reboot" -like de novo pattern formation 
from random initial conditions. 

Let us consider the simplified problem of regulating 
one domain, say the foot region versus the rest of the 
body. We consider this as a one dimensional problem 
as suggested by the well-defined head-foot-axis in Hydra. 
Fig. 1 formulates the target pattern of this problem. A 
basic means of communication between cells is local infor- 
mation transfer between neighbor cells through chemical 
agents, combined with information processing of a cell- 
internal gene network. Let us include local information 
transfer in the model, similar to the observed mechanism 
of direct contact induction in animals 0, . An in- 

teresting question is whether global positional informa- 
tion could be an outcome of local information transfer, 
in particular, when no concentration-dependent gene ac- 
tivation based on macroscopic gradients is involved. As 
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FIG. 1: The here considered morphogenetic problem of 
boundary formation and scale-independent proportion regu- 
lation of adjacent domains. The target pattern consists of 
one domain with a fraction of aN cells in state Ui = 2, a 
boundary state &i = 1 at the position aN, and cr; = for the 
remaining domain. 
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TABLE I: Rule table of the 3-state cellular automaton. First 
column: rule table index R. Second column: input states a = 
((Ti-i,(Ti,<Ti+i) at time t — 1. Third column: corresponding 
output states at time t. 



only further requirement, the up-down-symmetry of the 
body axis has to be broken locally in that case, e.g., by 
a spatially asymmetric receptor distribution [L4| . 

To define a model system that performs the above task, 
consider a one-dimensional stochastic cellular automaton 
[lH with parallel update. N cells are arranged on a one- 
dimensional grid, and each cell is labeled uniquely with 
an index i 6 {0, 1, N — 1}. Each cell can take n pos- 
sible states Ui <E {0, 1, .., n}. The state <7j(t) of cell i is a 
function of its own state <Ji{t — 1) and of its nearest neigh- 
bor's states <7£_i(i — 1) and <7j+i(t — 1) at time t — 1, i.e. 

<7i(t) = /(tTi_l(< - 1), CTi(t - 1), - 1)) (1) 

with / : {0,1, ...,n} 3 1— > {0,1, ...,n} (cellular automaton 
with neighborhood 3). At the system boundaries, for sim- 
plicity, we choose a discrete analogue of zero flux bound- 
ary conditions, i.e. we set <7_i = cat = const. = 0. Other 
choices, e.g. asymmetric boundaries with cell update de- 
pending only on the inner neighbor cell, lead to similar 
results. 

In the following, let us study the morphogenetic prob- 
lem formulated above for a — 0.3. We searched for pos- 
sible solutions in rule space by the aid of genetic algo- 
rithms (for details see [iff). The rule table of the best 
solution found is shown in Table I. Fig. 2 demonstrates 
the dynamics. Starting from a random initial configura- 
tion, the pattern self-organizes towards the target pattern 
within a finite number of updates. The finite size scaling 
of the self-organized relative domain size a as a function 




FIG. 2: A typical dynamical run for the cellular automata 
model, time runs from top to bottom. Starting from a random 
initial configuration, cells reorganize into an ordered expres- 
sion pattern corresponding to two asymmetric domains and a 
sharp domain boundary (deterministic dynamics, system size 
N = 150, color code of cell states: black a = 0, red at = 1, 
blue (Tj = 2) . 



of the number of cells N is shown in Fig. 3. In the limit 
of large system size, the boundary becomes sharp with a 
converging towards Ooo = 0.281 ± 0.001. The variance of 
a vanishes with a power of N, i.e. the relative size of fluc- 
tuations induced by different initial conditions becomes 
arbitrarily small with increasing system size. The pat- 
tern self-organization in this system is, therefore, quite 
robust. The main mechanism leading to stabilization at 
= 0.281 is a modulation of the traveling velocity of 
the right phase boundary in Fig. 2 due to particle inter- 
actions. On average, the boundary moves ca. one cell 
to the left per update step, whereas the left boundary 
moves one cell to the right every third update step. An 
intuitive concept for the dynamics of cellular automata 
phase boundaries views boundaries as moving particles. 
This so-called "particle computation" describes the phe- 
nomenology of complex cellular automata in terms of 
these soliton-like excitations [l7j . 

Let us now study the dynamics of the system under 
noise. For this purpose, we introduce stochastic update 
errors with probability p per cell, leading to an aver- 
age error rate r e = pN. In numerical simulations, we 
now have to apply a statistical method to measure the 
boundary position in order to get conclusive results also 
for high p: starting at i = 0, we move a measuring frame 
of w cells to the right and measure the fraction z of cells 
with eTj = 2 within the frame. The algorithm stops at 
some i when z drops below 1/2 and the boundary po- 
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FIG. 3: Finite size scaling of relative domain size a in equi- 
librium. In the limit of large system size, a converges to a 
fixed value Ooo = 0.281 ±0.001 (corresponding to line fit). In- 
set: finite size scaling of the variance Var[a(N)], fluctuations 
vanish with a power —1.3 of the system size. 




FIG. 4: Left panel: Cellular automata dynamics under noise 
(error rate r e = 0.005). Stochastic errors produce two differ- 
ent kinds of quasi-particle excitations, see right panels. Top: 
Interaction of the F-particle with the boundary readjusts the 
boundary two cells to the left. Bottom: A-particle moves 
boundary one cell to the right. 



sition is denned as i + w/2. It is easy to see that, for 
not too high p, there are only two different "particles" 
(i.e. state perturbations moving through the homoge- 
neous phases), as shown in Fig. 4. In the following, these 
particles are called V and A. The V particle is started 
in the a 2 phase by a stochastic error er^ = 2 — > =/= 2 
at some i < aN, moves to the right and, when reaching 
the domain boundary, readjusts it two cells to the left 
of its original position. The A particle is started in the 
(To phase by stochastic errors <7; = — > <7; 7^ at some 
% > aN and moves to the left. Interaction with the in- 
terface boundary readjusts it one cell to the right. Thus 
we find that the average position a* of the boundary is 
given by the rate equation 



2a*r e = (1 - a*) 



(2) 



i.e. a* = 1/3. Interestingly, for not too high error rates 
r e , a* is independent from r e and thus from p. Eqn. 
(2) also implies that the system undergoes a first order 
phase transition with respect to a* at p = 0; numerical 
evidence for this conclusion is given in [lfjj ] . Fluctuations 
of a around a* are Gaussian distributed with variance 
vanishing ~ iV _1 ^(|. The solution a* — 1/3 is stable 
only for < r e < 1/2. As shown in Fig. 4, the interaction 
of a r particle with the boundary needs only one update 
time step, whereas the boundary readjustment following 
a A particle interaction takes three update time steps. 
Therefore, the term on the right hand side of Eqn. (2), 
the flow rate of A particles at the boundary will saturate 
at 1/3 for large r e , leading to 2cv*r e = 1/3 with the 
solution 



O(N) 



(3) 



for r e > 1/2. Hence, there is a crossover from the solution 
a* = 1/3 to another solution vanishing with r" 1 around 
r e = 1/2 (Fig. 5). The finite size scaling term Q(N) 
can be estimated from the following consideration: for 
p — ► 1, the average domain size created by "pure chance" 

is given by a* = N^ 1 E^LoC 1 / 3 )" • n « ( 3 / 4 ) If 
the measuring window has size w, we obtain Q(N) ~ 
(3/4) w N- 1 . 

In a biological organism, a pattern has to be robust 
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FIG. 5: Average boundary position a* as an order parameter 
in the presence of noise. Transitions are shown as a function 
of the error rate r e . Numerics averaged over 200 initial con- 
ditions with 2 ■ 10 6 updates each. Dashed curves: mean field 
approximation given by Eqn. (3). Horizontal dashed lines: 
unperturbed solution a* = 0.281 and solution under noise 
a* = 1/3. 
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FIG. 6: Average boundary position a* as a function of error 
rate r e for different cell flow rates r/. Dashed curves: corre- 
sponding solutions of Eqn. (4). Note that a minimum error 
rate is necessary to regulate the boundary in the presence of 
cell flow. 

not only with respect to dynamical noise, but also with 
respect to "mechanical" perturbations. In Hydra, for ex- 
ample, there is a steady flow of cells directed towards the 
animal's head and foot, due to continued proliferation 
of stem cells. The stationary pattern of gene activity is 
maintained is spite of the cell flow. Let us now study 
the robustness of the model with respect to this type of 
perturbation. Consider a constant cell flow with rate r/ 
directed towards the left system boundary. In Eqn. (2), 
we now get an additional drift term ry on the left hand 
side: 2a*r e + rf = (1 — o*)r e , with the solution 

a * = { \ (l - £) if r e > r f and r e < 1/2 (4) 
\ if r e < r f 

with a* exhibiting a second order phase transition at the 
critical value r^ nt = rf. Below r^ nt , the domain size 



a* vanishes, and above r^" 11 it grows until it reaches the 
value a* max = 1/3 of the system without cell flow. The 
second order phase transition at rf lt bears some similar- 
ity with error catastrophes in models of viral evolution 
19]. Fig. 6 compares the numerical results with the mean 
field approximation of Eqn. (4). In numerical simula- 
tions, cell flow was realized by application of the transla- 
tion operator 8 (jj := Ui+\ to all cells with < i < N — 1 
every rj 1 time steps and leaving &n-i unchanged. Note 
that stochastic errors in dynamical updates for r e > rf 
indeed stabilize the global pattern against the mechanical 
stress of directed cell flow. 

To summarize, we considered a problem of pattern for- 
mation motivated by animal morphogenesis in a non- 
traditional setting. Accurate regulation of position in- 
formation, exhibiting proportional scaling with system 
size, and robust de novo pattern formation from random 
conditions have been obtained with a mechanism based 
on local information transfer rather than the Turing in- 
stability. Non-local information is transmitted through 
soliton-like quasi-particles instead of long-range gradi- 
ents, and fine-tuning of parameters is not needed. Noise 
contributes to the stability by generating quasi-particles 
that control the pattern. We observe considerable sta- 
bility also under cell flow. A first order phase transition 
is observed for vanishing noise and a second order phase 
transition at increasing cell flow. The pattern formation 
mechanism studied here is very general and not limited 
to cellular automata. In particular, implementations as 
regulatory networks work as well |l(j , and do not differ in 
complexity from regulatory circuits observed in the cell 
[Tl) . With this work we hope to inspire new approaches 
to biological pattern formation and to properties of non- 
equilibrium systems. 
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